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Abstract 

We present a perturbation method for determining the moment sta- 
bility of hnear ordinary differential equations with parametric forcing by 
colored noise. In particular, the forcing arises from passing white noise 
through an 71th order filter. We carry out a perturbation analysis based 
on a small parameter e that gives the amplitude of the forcing. Our per- 
turbation analysis is based on a ladder operator approach to the vector 
Ornstein-Uhlenbeck process. We can carry out our perturbation expan- 
sion to any order in e, for a large class linear filters, and for quite arbitrary 
linear systems. As an example we apply our results to the stochastically 
forced Mathieu equation. 

Subject Class: Primary: 93E15; Secondary: 60H10, 34D10. 

1 Introduction 

1.1 A Class of Stochastically Forced Linear Equations 

The original goal of this work was to develop a framework for analyzing the 
stability of the stochastically forced Mathieu equation: 

i + 7i+(w2^£/(t))x = 0, (1) 

where / is a stochastic process, and the stability is determined by the bound- 
edness of the second moment {{x'^{t))) [5l[T3]. Here, ((•)) denotes the sample- 
average. We wanted to avoid heuristic methods, and consider cases where f{t) 
is a stochastic process with a realistic power spectral density. In particular. 
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we do not want to assume that / is white noise. Hence we want to analyze 
the case where f{t) is colored noise. However, in order to rigorously derive 
a Fokkcr-Planck equation for a stochastic differential equation, the governing 
equation must include only white noise [S]. We can achieve both goals of rigor 
and realistic power spectral density by letting / be the output of a linear filter 
that is forced by a vector white noise ^. That is, 

s = Hs + ^(t), (2) 
/W = (a,s(t)), (3) 

where H is an nxn real, diagonalizable matrix, whose eigenvalues have negative 
real parts, a S M", and (•,•) is the standard inner product on C". Wc will 
take deterministic initial condition s(0) = 0. We assume the noise vector ^ is 
weighted white noise, meaning 

((a=0, Uit + r)em^B5{r), (4) 

where B is symmetric and positive semi-dcfinitc. Thus, when s(t) solves ([2]) it 
is a standard vector-valued Ornstein-Uhlcnbcck process. We refer to the scalar 
process, f{t) = (a, s(t)), as colored noise or as an nth-order filter provided 
s(t) solves Wc will make only mild requirements on the matrices H and 
B, thereby allowing for wide variability in the power spectral density of the 
resulting process (a, s{t)). Thus, in allowing for a wide range of choices of H, B, 
and a, our approach accommodates a broad class of colored noise forcing terms. 

In this paper wc will be concerned with the more general problem of linear 
equations that arc being parametrically forced by the function f{t) in equation 
That is equations of the form 

X = Aqx + £(a, s)AiX, (5) 

where s is the solution to the stochastic equation ([2]), x(t) G for some > 1, 
and Aq, Ai are N x N constant matrices. 

The purpose of this paper is to present a perturbation method (assuming e 
is small) for determining the stability of the solution x(t) of (O, by which we 
mean the boundedness of the second moments of x(t). However, our method 
applies to the pth moment, so we will not limit our analysis to second moments 
only. Van Kampen has presented a heuristic approach to the case of colored 
noise forcing, |24]. Though derived by completely different means, his result for 
the Mathieu equation ([1]) is the same as ours when considering only the first- 
moments, and without damping. He arrives at his result by truncating some 
series at order e^, and is only expected to be valid to this order. Our method is 
rigorous, can be applied to find solutions to any order in e, and applies to any 
moment. We discuss this further in §5.21 

We were originally interested in equation (IT|) as a model for the response of 
capillary gravity waves to a time-varying gravitational field arising from ran- 
dom vertical motions of a container with a free surface (as in [19]). Here f{t) 
represents the random fluctuations in acceleration. Since the Fourier transform 
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of an acceleration should vanish at zero, along with its derivative, the power 
spectral density of a realistic process / should satisfy 5(0) = S"(0) = 0. For 
example, we can construct a two-dimensional filter using the system ^ that 
has the power spectral density 



^(^) = L.2 , ..2^,^.2 , .,2V 



by choosing 



The formula for S'(w) in equation (jB]) follows from Corollary [7] in Appendix C. 

The stochastically forced Mathicu equation has been analyzed before, for 
instance in [31 [31 IH HI HH [12 but not for the case ([T]), or the general setting ([S]). 
In [ini [12 they consider additive forcing, and in [U [31 [1] they consider a different 
type of parametric forcing. In [9] they consider a different class of colored noises 
and study stability by truncating an infinite hierarchy of moment equations. 
Other studies concern Lyapunov stability, or rely on numerical methods. Our 
analysis applies to a broad class of equations ([5]) with a wide variety of forcing 
terms ([1]), is semi- analytical (relying only on numerics for the computation of 
eigenvalues of small matrices), and can be applied to any moment. 



1.2 Ladder Operators and the Vector Ornstein-Uhlenbeck 
Process 

Our perturbation analysis of the moment stability of equation ([S|) relies heav- 
ily on a simple characterization of the eigenvalues and eigenfunctions of the 
Fokker-Planck equation associated with equation ([2|). In particular, in ^and 
^we characterize the spectrum using ladder operators by generalizing Dirac's 
creation and annihilation operator approach to the quantum harmonic oscillator 

. An understanding of the spectrum and eigenfunctions in terms of its ladder 
operators is crucial to developing the perturbation theory in 21 Though other 
authors have used ladder operators for Ornstein-Uhlenbeck processes, they have 
only considered the scalar case n = 1 [2Ql El [23l [25] . We believe the extension 
to the vector case is by no means trivial, and is interesting in its own right. 

The probability density function P(si, . . . , s„, associated with the process 
s{t) defined by equations ([U and ^ satisfies the Fokker-Planck equation 

dtP = VP, with VP = idiv (BVP) - div (HsP) . (7) 

V is called the Fokker-Planck operator associated to ([U . See [H] for a deriva- 
tion of this equation. We note that the Fokker-Planck equation ([7]) is the same 
in both the Ito and Stratonovich interpretations because the matrix B is inde- 
pendent of s (see [H])- The operator V will play a crucial role in our stability 
analysis. 
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In ^J5]we begin by analyzing the operator D in terms of its associated ladder 
operators. That is, operators C satisfying the commutator equation 

[2?,/:]=/i£ (8) 

As in Dirac's theory of the harmonic oscillator, the significance of the ladder 
operators stems from the fact that if is an eigenfunction of V with eigenvalue 
A, then the function C(f> will either vanish, or be an eigenfunction of V with 
eigenvalue A + ^. 

In fJ5]wc show that we can construct the ladder operators by solving a matrix 
eigenvalue problem 

Ty = /iy, T = DA, (9) 

where A is an antisymmetric matrix and D is a symmetric matrix, expressed 
in terms of H, B. We show there are n raising operators, Ck, which generate 
new eigenfunctions of T) with an increase in the real part of the eigenvalue, 
and n lowering operators £-k that correspondingly decrease the real part of 
the eigenvalue. We also show that V can be expressed in terms of its ladder 
operators. In particular, 

n 

= (10) 

k=l 

where is the increment of the ladder operator Ck- That is, [D,Ck] = fJ-k^k- 
This representation is useful for determining the spectrum of T>. 
In ij3]we characterize the solutions of 

2?0 = x'/' (11) 

in terms of the ladder operators, C, and increments, fi, solving ([8|). In particular, 
we show that any eigenvalue x of 2? can be written as 

n 

Xk = - X! (12) 

where are the increments of the ladder operators with positive real parts, 
and the kj are non-negative integers. We will see that the increments fij are 
the negative of the eigenvalues of the matrix H defining the filter in equation 
([3]). We also show that any eigenfunction of V can be obtained by applying the 
ladder operator to the eigenfunction 4>o(s) associated with the eigenvalue x = 
of r>, which is the eigenvalue with the largest real part. 

The results summarized in the last paragraph rely on the fact that real parts 
of the eigenvalues of V are bounded above (see Lemma ^ , which is proved in 
pn [TSl [22], but we give a different and simple proof of this in Appendix B. 
Here, the domain of "D is the set of functions that have bounded moments of 
any order. The spectrum and eigenfunctions of T> have been studied before (see 
[T71 HSl [21] ) but not in the context of ladder operators. 
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1.3 Perturbation Expansion for Moment Stability Analy- 
sis 



In fj4] wc use the classical perturbation theory of eigenvalues to carry out an 
analysis of the stability of equation ([5]). Our analysis begins by considering 
the ODEs for s{t) and x(t) together as a single ODE system. The probability 
density function P{si, . . . , s„, xi, . . . , xjy, t) for the combined system ^ and (O 
solves the Fokker-Planck equation 

dtP = idiv, (BV,P) - div, (HsP) - div, ((Ao + e(a, s)Ai)xP) (13) 

The notation divg and Vg refer to divergence and gradient with respect to only 
the Sj variables, and similarly divx is divergence in x variables. Equation (jl3p 
is the same in both the Ito and Stratonovich interpretations because the matrix 
B is independent of s and x (see |12)). 

We can derive an equation for the pth marginal moments by multiplying (jl3p 
by monomials x" and integrating with respect to dx, where a is a multi-index 
of order p. The result is an equation for m(s,<), a vector of the pth marginal 
moments, which is of the form 

atm = 2?m + rom + e(a,s)rim. (14) 

Note that 2? is a differential operator in the s variables only, 

= idiv, (BWsV) - divs (Hs(/j) . (15) 

In equation ([T^ each component of m(s, t) is of the form Jj^^ x"P(x, s, t)dx for 
some multi- index a with |a| = p. 2?m indicates V applied to each component 
of m. For much of our analysis we can assume that the matrices Tg and Ti are 
given to us, but we illustrate how to obtain these matrices from the matrices 
Ao and Ai for the particular case of the Mathieu equation in ^ The matrices 
To, Ti in ()14p are constant and depend on which moments one is considering 

(see example in equation (j53p ). There are J = ^ ^ ^p ^ ^ distinct pth 

order monomials in N variables, therefore To and Fi are J x J matrices. 

As in a standard stability analysis, in order to determine the stability of 
P^ . we look for solutions of the form m(s,i) = e'^*m(s). Our equation for 
m(s) becomes 

Am = 2?m + rom + £(a,s)rim. (16) 

That is, the equation for the pth marginal moments of x{t) can be written as 
an eigenvalue problem, and stability is decided by the sign of the real part of 
the largest eigenvalue. 

We do a perturbation analysis assuming that the magnitude e of the forcing 
is small. Our analysis relies on the fact that that when e = the eigenfunc- 
tions of equation (fT6)) are the direct product of the eigenfunctions of V and the 
eigenvectors of the matrix Fq. 
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A key observation (see Lemma [H]) for the perturbation analysis is that for 
any vector a we can determine constants ak and /3k such that (a, s) can be 
written as 

n 

(a,s) =^K/:fe + /?fc/:_fc), (17) 

k=l 

where £±k are the ladder operators satisfying ([5]). The proof of Lemma [5] is 
given in Appendix D. 

In 21 we show that when e = 0, the eigenvalue of equation ([T5| with the 
largest real part is the same as the largest eigenvalue of the matrix Tq. If Aq is 
this unperturbed eigenvalue, then A(e) = Aq + A2£^ + . . . with 

J 

A2 = ^(V'i,ri</.^.)(V„ri0i)G(z/i -i.,) (18) 

where Vi = Aq, i^j are the eigenvalues of To, and 4>ji '4>j are the eigenvectors and 
normalized adjoint eigenvectors of Fq. Equation (jl8p uses the extended power 
spectral density G{z), which is defined for a general stationary random process 
in ((79| . and is given explicitly in (|8T|) for the filter (a, s{t)). 

The form of A2 in is derived for forcing terms that have the form 
however, the fact that this is simply a weighted sum of values of G, whose 
coefficients depend only on ro,ri (which do not depend on the filter), suggests 
such a formula could hold for any process with a well-defined extended power 
spectral density. We have carried the perturbation analysis to higher orders, 
but the higher-order coefficients do not appear to have such a simple form as in 
equation (fTS]) . 

The method in f}4] involves constructing matrices Fq and Fi, which, as men- 
tioned earlier, depend on Aq, Ai, and the representation of the pth marginal 
moments as a vector. We use the stochastic Mathieu equation as a specific ex- 
ample in ^JS] In i j5.3l we discuss a numerical method for determining the stability 
of (O without assuming that e is small. We compare these numerical results 
to our perturbation results up to both second and fourth order for the Mathieu 
equation, and show they are in excellent agreement. In ij5.4l we give a second 
representation (whose derivation is given in |B]) for A2 that does not involve the 
matrices Fq and Fi, but deals directly with the matrices Aq and Ai. We have 
found that this representation simplifies numerical computations. 

2 Existence and Properties of Ladder Operators 

In this section we define the notion of a ladder operator, show how to construct 
these operators, and prove some basic lemmas about them. Lemma [1] shows 
how ladder operators can be used to generate new eigenfunctions that have their 
eigenvalue changed by the increment of the ladder operator. Lemma [5] shows 
how to find the ladder operators Ck and their increments by solving a matrix 
eigenvalue problem. Lemma[3]shows that the increments of the ladder operators 
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are zero and ±Hk ,where —^k are the eigenvalues of the matrix H defining our 
filter (see equation ([2]))- Lemma |4] gives the commutator relations for the ladder 
operators, and Lemma [S] shows that the operator 2? can be expressed as a 
weighted sum of C-kCk, where Ck are the ladder operators. Throughout this 
section (and the rest of the paper) the operator D is that defined in . 

The basic lemmas in this section are crucial to the rest of this paper, and 
hence we have tried to write this section so that the lemmas stand out clearly. 
Though the lemmas are all easily stated, the proofs of some of the lemmas are 
quite technical, especially when attention is given to ensuring that they apply 
for complex eigenvalues of the matrix T. For this reason, we have relegated 
many of the proofs to Appendix A. 

Before discussing ladder operators it should be noted that we define the 
domain of "D as the set of functions that have bounded moments of any order. 
Thus, our definition of the domain of V differs from that given in [T7]. In 
that paper they defined the domain based on the exponential decay of the 
eigenfunction $o that we define in equation ([5^ and discuss in 321 The two 
definitions of the domain give the identical eigenfunctions, but we believe ours 
is more natural since it does not require knowing the solution ahead of time. In 
[17] they discuss a continuous spectrum that arises if the domain is defined so 
that the eigenfunctions of V are only required to be square integrable (or some 
similarly less restrictive condition). An examination of these eigenfunctions 
shows that they have a power law decay as s goes to infinity, and hence do not 
have moments of all orders. Hence our definition of the domain also excludes 
this continuous spectrum. 

We now give a definition of a ladder operator of V. 

Definition 1. An operator £ is a ladder operator for V with increment fi if 
[D, C] = jiC for some e C, where [• , •] denotes the commutator [D, C] = 

vc - cv. 

The following lemma shows that ladder operators can be used to generate 
new eigenfunctions from ones that we already know. 

Lemma 1. Suppose C is a ladder operator such that [D^C] = fiC. Let (p he 
an eigenfunction of T) with eigenvalue x- Then either Ccj) = 0, or C(j) is an 
eigenfunction of V with eigenvalue x + A^- 

Proof. We have VCcj) — CD^) = ^Ccf). Since (j) is an eigenfunction of P, this gives 
us VCcj) ={x + □ 

We defined the domain of T) to be the set of functions that have moments of 
all orders. It should be noted that Lemma [T] would not apply if the domain had 
been (for example) the set of all square integrable functions. In that case a third 
possibility would exist. It could be that the function (p is square integrable, but 
the function C(j} is not. Thus C(j} would not generate a new eigenfunction. 

We will show that T) has 2n + 1 ladder operators £±fc, k = 0, . . . ,n. We 
begin by decomposing T) into simple differential and multiplicative operators. 
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Definition 2. We define the operators Lj,j = 1, . . . , 2n + 1 as follows. 



Lj(l) = ds^4> forj = l,...,n 

L.j+n(f) = Sj(j) for j = 1, . . . , n (19) 

L2n+i't> = I4> for j ^2n + l 

Here I is the identity operator. Note that [Lj,Lk\ = unless |j — fc| = n, 
and [Lj , = T. In particular, we have 

[ij, Lfc+„] = (5jfel j,k^l,...,n (20) 

We note that D can be expressed in the operators Lj as 

2n+l ^ 

^ = E 2'^JkL,L,. (21) 

We let D denote the symmetric matrix with components djk in ((2T|) . The choice 
of djfe in ()2ip is not unique, but we make an explicit choice that makes this ma- 
trix symmetric. If we let A denote the antisymmetric matrix with components 
ajk given by 

[L,,Lk]^a,kI, (22) 
then we have explicit expressions for A and D 

0„ I„ \ / B -H 

A= I I„ 0„ ; I D= o„ ; I- (23) 

••• / V ... -tr(H) 

For details regarding the construction of D see Lemma [HI in Appendix A. 

Just as 2? has a representation in terms of the operators Lj , its ladder oper- 
ators will also be expressed in terms of the Lj. Consider an operator 

2n+l 



We write y for the vector of coefficients of C. From the representations (|2ip 
and ([24|) . we see that the commutator [P, C\ involves sums of terms of the form 
LiLjLk — LkLiLj, which do not at first sight appear to be linear in the operators 
L„n m = 1, . . . , 2ri + 1. However, by twice applying the commutator relations 
in equation (|20p we can show that LiLjLk — LkLiLj is in fact a sum of the Lm- 
Determining the coefficient vector y and increment ^ thus becomes a matrix 
eigenvalue problem. The details of how we arrive at this form are given in 
Appendix A. Here we will merely state the result of these manipulations. 

Lemma 2. //y is the vector of coefficients for C, as defined as in equation {24% 
then the equation [D, C] ~ fiC can be written as a matrix eigenvalue problem 
Ty = /iy, where T = DA, and D and A are defined in ([23]). 
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We make the assumption that the eigenvalues of H have negative real parts, 
and the eigenvectors form a complete set. For simplicity of the arguments, we 
will also assume that the eigenvalues of H are simple. By explicitly writing out 
the eigenvalue problem Ty = /iy we can determine the eigenvalues in terms 
of the eigenvalues of the matrix H. We will give the details of the proof in 
Appendix A. 

Lemma 3. The eigenvalues ofT — DA are {0, ±/Xfc}, k — 1, . . . , n, where ~fik 
are the eigenvalues of the matrix H. 

Note that Cq is the identity operator with increment 0. Thus, our analysis 
only involves the 2n ladder operators C±k for fc = 1, . . . , n. 

In doing the perturbation expansion it will be necessary to have the com- 
mutator relations of the operators Ci. Finding the commutator relations for 
[Cj^Ck] can be turned into a linear algebra problem involving the eigenvectors 
of the matrix T. In particular, using equations and ([2^ we get 

[£„£fc] = (yjAyfc)X. (25) 
From equation ([^5)) it is easily seen that 

AD = -(DA)^ = -T"^ (26) 

If DAy = ny, then multiplying both sides of this by A and using equation 
([M]) we see that Ay is an eigenvector of T-^ with eigenvalue —fi. With this in 
mind, the right hand side of equation can be written as the inner product 
between the vector and the adjoint eigenvector of T associated with — /Xfc. 
Using the fact that the eigenvectors and adjoint eigenvectors of a matrix form 
a bi-orthogonal set, we can arrive at a simple expression for the commutators. 

When dealing with complex quantities, the notation in this argument gets 
to be a bit tedious, and we will leave the details to Appendix A. The final 
commutator result is given by the following lemma. 

Lemma 4. For j,k > I we have [Cj, Ck] ~ and [C-j, Ck] ~ Sjkl- 

In Dirac's theory of the harmonic oscillator, he shows that the Hamiltonian 
operator can be written as the product of the raising and lowering operators. 
We now generalize this result to the vector case. In this case the operator V 
can be written as a weighted sum of the products of the raising and lowering 
operators. The next lemma shows that the weights are in fact the eigenvalues 
/Life of the matrix T. We leave the proof of this lemma to Appendix A, but note 
that its proof is probably the most subtle one in this paper. 

Lemma 5. The differential operator T> ~ jii \'^i-3^i^j ^'^'"^ written as 

n 

V = J2t^kC-kCk. (27) 

fc=i 

An important feature of the decomposition p7|) is that only terms of the 
form fc > 0, appear (there are no terms of the form CkC-k for fc > 0). 
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3 Eigenvalues and Eigenfunctions of V 



In this section wc will use the ladder operator formalism to completely charac- 
terize the eigenvalues and eigenfunctions of the operator V. We note that the 
spectrum of V has already been studied and characterized [T71 [THl [H] , but not 
in terms of ladder operators. We include another proof of those results because 
the characterization in terms of ladder operators is used in the perturbation 
analysis in 21 

As with Dirac's theory of the quantum harmonic oscillator, the analysis of 
the spectrum using ladder operators requires that the real part of the spectrum 
is bounded above. We will now state this as a lemma, but leave the proof to 
Appendix B. 

Lemma 6. The real part of spectrum of the operator T), as defined in (jlSp . is 
bounded above. 

The following theorem will allow us to characterize the eigenfunction asso- 
ciated with the eigenvalue with the largest real part. 

Theorem 1. Let <&(s) be an eigenfunction ofD (as in equation (|15[) ) associated 
with the eigenvalue having the largest real part. We must have £/c<i? = for 
/c = 1, . . . , n. 

Proof. Suppose $(s) is an eigenfunction of T) with eigenvalue x- = Ck^ ^ 
0, then 4'(s) will be an eigenfunction of V with eigenvalue x + Mfe- This will 
give us an eigenvalue with a larger real part than x- Hence if x is the eigenvalue 
with the largest real part, then = for all k. □ 

Remark 1. The system £/j"I> = for each A: = 1, . . . , n is an over-determined 
system of first order differential equations. The fact that a solution exists is 
non-trivial. However, the fact that [Ck^Cj] = implies that the Frobenius 
Theorem applies (see [11 E]), which guarantees the system is solvable. 

Remark 2. As in the comment following Lemma [l] we should note that the 
domain of T) is defined as the set of functions that have moments of all orders. 
If the domain of T) were defined using the less stringent requirement that the 
eigenfunctions were square integrable, it would not be necessary that Ck^ ~ 
for all k. This is because in this case Ck^ does not have to generate a new 
eigenfunction. It could instead produce a function that is not square integrable. 

By Theorem [U the "top" eigenfunction $o(s) (i.e. the eigenfunction associ- 
ated to the largest eigenvalue of V) must satisfy the equations £^$0 = 0. If 
is the eigenvector of T associated with the eigenvalue fjLk, and if fik 7^ 0, then 
the last component of y'^ vanishes (see the proof of Lemma [3] in Appendix A) . 
That is, we can write 




(28) 
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Using equation 1^^, and the definition of the operators Lk in the equations 
-Cfc^o = can thus be written as 




fc = l,...,n (29) 
is-^Ss), then equations ([^^ wiU 



(30) 



where 



P = [pi,p2,...,p„], 



Q = [qi,q2,...,qri]. 



(31) 



If P is invertible, this gives us S = (P-^)~^Q-^. It is not clear that P is 
invertible, or that S is symmetric. However, under certain weak assumptions 
on H and B (see Definition |3] and Lemma [7] below) this will be the case. If these 
assumptions hold, it is convenient to write = PQ^^. We now define the 
notion of a controllable pair. 

Definition 3. The matrices H and B will be said to form a controllable pair 
if there is no nontrivial vector z such that z-^H'^B = for fc = 0, . . . , n — 1. 
This is equivalent to requiring rank C = n, where C is the n x matrix 



In Appendix B we prove the following lemma. 

Lemma 7. Assuming all of the eigenvalues o/H have real parts less than zero, 
the eigenvectors of H are complete, and that B is positive semidefinite, then 
= PQ~^ is symmetric and positive semi- definite. //H and B also form a 
controllable pair, then the matrix is positive definite, and hence the matrices 
S = QP^^ and P are non-singular. 

Requiring (H, B) to be a controllable pair eliminates some "degenerate" 
types of filters. For instance, if H = diag(— /ii, — ^2) and B = diag(l,0) 
then (H,B) is not a controllable pair. In this case, si{t) is a scalar Ornstein- 
Uhlenbeck process, but S2(t) is deterministic, so s{t) is not a genuine two- 
dimensional Ornstein-Uhlenbeck process, but rather it is a one-dimensional pro- 
cess with an appended deterministic component. 

Definition 4. We will say that the n x n real matrices H and B satisfy the 
basic conditions if 

(i) B is symmetric and positive semi-definite 

(ii) H has simple eigenvalues {~f^k}]^^i with Re [/i^] > for fc = 1, . . . , n 

(iii) (H, B) form a controllable pair (Def. [3]). 

The requirement of simple eigenvalues for H is for convenience and could be 
replaced with the requirement of a complete set of eigenvectors. 



C = [B, HB, . . H"-iB]. 
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Lemma 8. Assuming H and B satisfy the basic conditions (Def. the eigen- 
value Xo with the largest real part of V is simple, and the eigenfunction <I>o (s) 
associated with it is given by 

1>o(s)=exp(^-i(s,£s)^ , (32) 
where S = QP""'^. Moreover, xo ~ 0- 

Proof. Without loss of generahty we look for solutions of the form $o(s) = 
g-i/'o(s)_ order to satisfy equations ([^^ we must have 

Pfc • VVo + (qfe • s) = (33) 

A direct calculations shows that ipo = — ^(s, Ss) satisfies this equation. If we 
have another solution to this equation, say ipi, then the difference V'o — i^i 
between these solutions will satisfy • V('0o ~ V'l) — 0, for k = 1, . . . ,n. The 
vectors are complete (they are the eigenvectors of H^), which implies that 
ipQ — -01 is a constant. This in turn implies that the eigenfunctions associated 
with each of the solutions ipoi^) £ire multiples of each other, hence xo is simple. 

From Lemma [SI the real part of the the spectrum of V is bounded above. 
From Lemma [5l we can write V = l^-k^-k^k- Hence, when we apply V 

to the eigenfunction $o(s) associated with xo, we will get I?$o = because 
'Cfc^o = 0, fc = 1, . . . , n, by TheoremUl Hence = 0. □ 

Theorem 2. Let H and B satisfy the basic conditions (Def. x is an 

eigenvalue of V if and only if it can be written as in equation where kj 

are non-negative integers, and —^Lj are the eigenvalues ofH. The eigenfunction 
associated with Xk «s given by $k(s) = £'!_\£j^2---'^-n*J'o(s)- where $o(s) is 
defined as in equation i32\) . 

Proof. From Lemmas [5] and [51 the real part of the the spectrum of V is bounded 
above by 0, and x = is an eigenvalue of V that has the form (|12p . If x is any 
other eigenvalue, and $ is its eigenfunction, then there must be at least one value 
of k such that Ck^ ^ 0. If this new eigenvalue has the form given in equation 
(fT2)) . then the previous one will too. We can keep carrying out this process 
obtaining eigenvalues with larger real parts. This process must eventually end 
since the real part of the spectrum is bounded above. The only way it can end 
is when we arrive at the largest eigenvalue, which we have already seen, is zero. 
This implies equation 

The argument in the last paragraph shows that any eigenvalue of 2? must 
be of the form p2p . To show that any number Xk of this form must be an 
eigenvalue of V wc show that $k(s) = £?i\£^i^2---'^-"n'^o(s) is the eigenfimction 
associated with Xk- This follows from Lemma [T^ in Appendix B. □ 

4 Perturbation Method 

The marginal-moment equation (|14p is derived by multiplying (jl3p by a mono- 
mial x" for some multi-index a, then integrating with respect to x. If this is 
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done for each multi-index of order p, we derive a set of equations for the pth 
marginal moments. If we collect the pth marginal moments into a vector m, 
we arrive at ([T?]). The matrices Fq, Fi depend not only on Aq, Ai, but also on 
our mapping of the pth marginal moments into m. For this reason, we do not 
write the explicit form of Fo,Fi in this section, but we do write them out for 
the example of second marginal moments for the Mathieu equation in |J5] 

We let cf)j denote the eigenvectors of Fq with eigenvalues ja,- . Wc let be 

T 

the normalized adjoint eigenvectors, so that tpkCpj = (i/'j., </)j) = Skj- We may 
assume without loss of generality that the are ordered so that Re [i^i] > Re [vj] 
for all j. 

We expand the unknowns as scries in e, 

A = Ao + eAi + £^A2 + . . . , m(s) = mo(s) + emi(s) + e^m2(s) + . . . , (34) 

and solve for the terms of these series. If we substitute these expansions into 
(fTB|) . and collect the zeroth-order terms, we get 

Aomo = Vmo + Fomo. (35) 

The eigenfunctions of V are scalar- valued, and the eigenvectors of Fq are con- 
stant vectors. Assuming that both the eigenfunctions of V and the eigenvectors 
of Fq are complete, then the most general solution mo to ((35|) will be a product 
of an eigenfunction of I? with an eigenvector of Fo , and Ao will be the sum of 
the eigenvalues of D and Fq. We are interested in the largest eigenvalue, so we 
take 

mo(s) = $o(s)0i, (36) 

and Ao = I'l because is the largest eigenvalue of V and V^q = (Lemma [8]) , 
and vi was selected to have the largest possible real part (note that the choice 
of vi need not be unique). 

The form of the forcing in (O allows us to represent (a, s) in terms of the lad- 
der operators. In particular, the parametric forcing by the linear filter (a, s(f)) 
results in the presence of the first-order polynomial (a, s) in the Fokker-Planck 
equation, and thus to the term e(a, s)Fim in the moment equation (jl6p . Since 
the ladder operators, C±k, are linear combinations of first-order operators 
and monomials Sj , it is reasonable to try to write (a, s) as a linear combination 
of C±k- The completeness of the eigenvectors of H allows us to do this, greatly 
simplifying our perturbation analysis. 

Lemma 9. If the eigenvectors o/H are complete, and a^, h o,'"'^ defined as in 
equation ([55)) (Appendix D), then 

n 

M^Y.^akCk+hC-k). (37) 
fe=i 

The proof of Lemma [9] is in Appendix D. Formula (|37p ensures that the 
coefficients and will appear in the coefficients of the perturbation expan- 
sions ([M]) . Wc show in Appendix C that the extended power spectral density, 
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G{z), of (a, s{t)) can also be expressed in terms of ak and /3fc (TheoremlHl). This 
allows us to derive a simple formula for the order coefficient of A(e) in terms 
oiG{z) (Theorem O. 

f N + p - 1 \ . . 
Recall that there are J = distinct nth order monomials in 

\ P J 

N variables, and that Tq and Fi arc J x J matrices. We will assume that Tq 
has a complete set of eigenvectors, which is the case for the Mathieu equation, 
and occurs whenever the eigenvectors of Ag are complete. The following lemma 
gives solvability conditions that will be used repeatedly in our analysis. 

Lemma 10. Let H and B satisfy the basic conditions (Def. Suppose that 
Tq has a complete set of eigenvectors {4>j}j^i, with eigenvalues Vj, normal- 
ized adjoint eigenvectors {'0j}/=i; that $(§) is an eigenfunction of D with 
eigenvalue —p.. If n ^ 0, then then the equation 

(Ao-2?-ro)m = $(s)b 

has a solution given by 

m(s) = <l>(s)^ ; 0^. (38) 

On the other hand, if fi = (and hence $(s) = $o(s)j and vi ^ Vj for j > 1, 
then the equation 

(Ao-2?-ro)m = $o(s)b 
has a solution if and only if {xpi,h) =0. In this case, the solution is 

m(s) = K$o(s)</.i + $o(s)V^^^^^0j. (39) 

where k is an arbitrary constant. 

The constant k can be used to choose a normalization for m. Wc do not 
need to choose a specific normalization for m, so we set k = because it is 
convenient. One can check that if Aq has a complete set of eigenvectors then 
Tq will too. 

Proof. If /i 7^ 0, then when we write b in the (pj basis, and make the ansatz 
m(s) = $(s)c, where c is a constant vector, we arrive at the expression for 
m in equation (|38)) . If /x = 0, and hence ^(s) = $o(s), then we cannot solve 
this equation if b has any component in the direction of </)j. This gives the 
compatibility condition {xpi,h) ~ 0. Assuming this holds, the solution is given 
by equation ([39|) . □ 

We will now describe the outline of the perturbation analysis. In order to 
help us describe the perturbation analysis we will use the following definition. 
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Definition 5. We say a function f(s) is in Vk if it can be written as the sum 
of eigenfunctions of T> times constant vectors, where each of the eigenfunctions 
is the product of k or fewer ladder operators — 1, . . . ,n applied to the 

eigenfunction $o(s)- 

The following lemma will be used in our perturbation analysis. 

Lemma 11. If h{s) G Vk, then g{s) = (a, s)/i(s) is in Vk+i- 

Proof. This is almost a direct consequence of Lemmas |4] and [9] From Lemma 
[9] we know that g(s) can be written as a sum of terms involving C-jh{s) and 
Cjh{s) where j > 0. By definition, each of the terms C-jh{s) are in Vfe+i. On 
the other hand, the commutator relations [Cj,C-k] = ~SjkI from Lemma |31 
and the fact that £j^o{s) = for j > 0, can be used to show that Cjh{s) is in 
Vk-i That is, Cj has either canceled out a previous term C-j applied to $0, or 
it commutes with all of the previous operators applied to $o- yielding the zero 
function because Cj^o = for j > 0. □ 

The perturbation analysis proceeds as follows. We have a zeroth-order so- 
lution mo = <Pi^o, which is clearly in Vq. We will see by induction, that the 
function mfc(s) will be in Vk- 

The equation at each higher order will be of the form 

(Ao - I? - To)mk = -Afemo + r^^i (s) (40) 

where rfc_i(s) is function that can be computed using the m^ and Xj for j < k. 
In particular, we have 

fc-i 

rfe_i(s) = - '^Xjirik-j + (a,s)rimfc_i 

j=i 

Assuming that for j < k the functions mj(s) are in Vj, then Lemma [Til ensures 
that the term rfc_i(s) will be in Vk- We can write 

rfe_i(s) = $o(s)bfc_i +ffc_i(s) 

where the term ffc_i(s) can be written as a sum of eigenfunctions of V times 
constant vectors, where none of the eigenfunctions is $o(s). With this in mind 
we use Lemma [TU] to see that we will be able to solve equation if and only 
if \k{'4>i,(f>i) = (■»/'i,bfc_i), and hence 

Afe (■0i,bfc-l). 

Once we have chosen A^ in this way, we can solve for m^, and it will clearly 
be in V^, thus allowing us to continue the process to the next value of k by 
induction. 

Terms in rfc_i(s) proportional to $0 can only arise at even steps in the pro- 
cess (i.e. equations for A2j, m2j) because £fc$fc = — $0 (see equation |42)) . These 
terms proportional to $0 must satisfy the compatibility condition {xpi,h) = 
as in Lemma ITOl 
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4.1 First Order 

To simplify notation, we make the following definition. 

Definition 6. The funetions $fe(s) are defined as 

$fc(s) = /:_fc$o(s), k = l,...,n, (41) 

where ^o{s) is the eigenfunction of T) associated with the eigenvalue with the 
largest real part. Note that, from Lemma H] and Theorem [U we have 

£fc$fc(s) = -$o(s), fc-l,...,n, (42) 

because = AX-fe^o = {C-kCk - l)*o = -$o- 

Substituting p4p into ([T6|) and collecting terms of order e, we get the equa- 
tion for mi 

{Xo-V- ro)mi = -Aimo + (a, s)rimo. (43) 

It is not hard to show that the eigenvalue A(e) = Aq + eAi + . . . must be an even 
function of e. This is also intuitive because the sign of e plays no role in ([S]). 
Thus, it is no surprise that Ai = 0. 

Lemma 12. //H and B satisfy the basic conditions (Def. we have Ai = 
and 

n 

mi = ^$fc(s)cfc (44) 

fc=i 

where <^k is defined in (|4T|) . and 

Ck = 2^ cj) (45 

Proof. Using ^7}i and that Ck^o = 0, £-/£$o = $fe, we have 

n 

(a,s)rimo = ^/3fcri</.i$fe(s) (46) 

k=l 

Thus, the right side of P5)) is a finite sum of terms proportional to ^o,^i, . . . ,^n, 
and each term can be treated separately. We now apply Lemma [TU] to equation 
using equation (|46p . The only term proportional to $o is — Aimo. But ac- 
cording Lemma [TO] this means {ipi, Xicf)^) = 0. Hence, Ai = 0. The expression 
in (p8| applied to the terms for fc > gives the expression in (|44)) . □ 

4.2 Second Order 

Substituting ([M)) into (|16p and collecting terms of order e^, we get the equation 
for m2 

(Ao - X> - ro)m2 = -Aamo + (a, s)rimi. (47) 

The situation here is similar to that for mi , except that the terms proportional 
to $o(s) come from mo as well as terms of the form £fc$fe(s) = — $o(s). 
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Lemma 13. //H and B satisfy the basic conditions (Def. the compatibility 
condition for ( |^7| j implies 

n 

X2 = -Y,{ipi,akriCk}, (48) 

fe=i 

where Ck is defined as in 

Proof. Lemma [T^ shows that mi = X]fe=i ^fc(s)cfe. This fact, and an applica- 
tion of the result in Lemma [5] implies that 

(n \ n / ^ \ 

J2 (^i^i + Pi^^i Ti $fc(s)cfc = -$o(s) "fcriCfc +. . . 
1=1 / k=l \k=l / 

where the term on the right is the only term proportional to $o- We used (|42p 
to write £fc$fc(s) = CkC-k'^ois) = -$o(s). 

Using the form of mo in ((36)) . the compatibility condition from Lemma [TOl 
implies -A2(V'i,0i> = i'^i^Y.k^i'^k^i^k) ■ Hence 

n n 

A2 = -(V'lj^afcriCfc) = - ^(V'i,afcriCfc). 
fc=l fc=l 

□ 

Computing the expression for m2 is a simple exercise, but we do not write it 
here. Continuing this process for higher order terms is straightforward, though 
grows more tedious with each successive order. 

Lemma [T2] allows us to compute A2 , but a nice feature of the second order 
term A2 , is that it can be expressed by a simple formula involving the extended 
power spectral density G of the process (a, s(t)) (see Appendix C). We prove 
the following theorem in Appendix D. 

Theorem 3. // Re [vi — Vj + /.jfe] > for each j ~ 1, . . . , J and k ~ 1, . . . ,n, 
and i/ H and B satisfy the basic conditions (Def. then 

.7 

A2 - 5^(V'i,ri</,^.)(t/>^,ri(/.i)G(z.i - V,). (49) 

Here G{z) is the extended power spectral density of the forcing term (a, s) . 

Remark 3. Note that the coefficients {'ipi,Ti<pj){ipj,Ti<pi) and the differences 
vi — Vj depend only on the differential equation for x (i.e. only on the matrices 
Aq and Ai), and the function G depends only on the filter (a, s(t)) (i.e. on 
H, B, a). It would be interesting to investigate whether the same form as in (|49p 
would hold for any asymptotically stationary filter. That is, if the expression for 
A2 would be a linear combination of values of G, where the coefficients depend 
only on the physical system, and the places where G is evaluated are given by 
the eigenvalues of that system. 
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5 Applications 

5.1 Second Moments for the Mathieu Equation 

We can write the Mathieu equation ([1]) as in ((Sj using a two-dimensional vector 
X"^ {xi,X2). In this case the matrices Aq, Ai in equation ([5]) are 



We will consider the stability for the second moments. We define 

TOjfe(s,i)= / XjXkP{xi,X2,s,t)dxidx2, (51) 
In this case The Fokker Planck equation (fT3|) can be written as 



dtP ^VP- — ix2P) - ^ {{-ujIxi - JX2 + (a, s)a;i) P) (52) 

If we multiply equation (j52p by x^, and integrate over all values of Xi and 
X2, after integrating by parts we get the equation 

dtmii ~ Dmii + 2toi2 

Similarly multiplying equation (|52p by 2:1X2 and 2:2, integrating over all xi and 
0:2, and applying integration by parts, we get the equations 

dtmi2 = 2?TOi2 - ^Imii - 7777.12 + m22 + (a, s)mii 

and 

dtm22 = 'Dm22 - 2wQmi2 - 2777122 + 2(a, s)mi2. 

If we let m = (77711,77112,77122)"^ this can be written in the form of equation 
where the matrices ro,ri in are given by 



(53) 



After assuming temporal behavior of the form e^^ we arrive at the eigenvalue 
problem 

Am = Pm + Tom + e(a, s)rim. (54) 

for m(s) and A, which is the same as ([TS]). We will now apply the results of 
Theorem |3] to this set of equations. 

In the case of second moments, the eigenvalues Vj of Tq are given by sums 
of two eigenvalues of Aq. I.e., Vj = cr^ + Om where Uk are eigenvalues of Aq. 
In the case of the Mathieu equation, the eigenvalues of Aq are cri, (J2, where 

CTi = CT2 = ^ ^ Hence, there are three choices of t^i, given by vi ~ —7, 






2 






r 

























-2^0 


-27 / 




^0 


2 
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or j/1 = — 7 ± i\/ Abj^ — 7^, since they all have the same real part. In the case 
v\ = —7, we have (t/j^jFi^^) = 0, so the G(0) term does not appear. We also 
have (V'i,ri02)(i/'2,ri0i) = (i/)i,ri03)(V'3,ri0i) = s^^- Hence 

(G(i^i - v^) + G(z.i - = ^^^1^2 ^ (^^^40.2 - 72^ , (55) 



4wo - r 

where S'(w) is the power spectral density of (a, s(t)). This follows because 
(without loss of generality, taking vi = — 7— i-y/4a;g — 7^ = 173) we have z/i — 1^2 = 
z ^ 4wq — 7^ = —{v\ — t^a), and G{iuj) + G{—iLu) = ^(u;) (see Appendix C). 
If we take either i/i ~ —7 ± ^ Asjj^ — 7^, then the expressions for A2 are 

2 



4cj2 - 72 
2 

4cj^ - 7^ 



G (^z\/4c.2-72^ _ 2G(0) 



(+) A2 = 



(-) ^2 = x3^(G(-*V4a.o'-y )-2G(0) 



Both cases have the same real part of A2 



1 



[^2] = -2—1 ( S ( V4^§ - 7^ ) - 2^(0) 



4w^-72 V / 7 

which is less than the expression in (j55|) . Hence, we have proved 



Theorem 4. //H and B satisfy the basic conditions (Def. then the second 
moments of the Mathieu equation flp become unstable when A(e) > where 



5.2 Comparing Moments for the Mathieu Equation 

If we perform the same analysis as in ijS.li but for the first and third marginal 
moment equations instead of the second marginal moment equation, we obtain 
results similar to Theorem S) If we denote the largest eigenvalue of the pth 
moment operator 2? + To + e{a, s)Ti as X^p\ then up to second order, we have 



AW(£) 

A(3)(£) 



-7 + *v/4^o'-72 , /G(*V4c.2_^2) G{0) 



Alu^ - Y 4iu5 - 72 



£2 



-37 + ^40;^ -7^ 



3G 



(-^v/4^g-7^) 4G(^V4a.g-7^) g(o) 
4w^ -72 ^ 4a;^ - 7^ - 7^ 
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(Fi and Fq depend on p, but we do not make that explicit in our notation.) It 
is only the real parts of the eigenvalues that factor into the stability. We have 



Re 
Re 
Re 



In [21], there is a heuristic treatment of the first moments of x(i). There, Van 
Kampen writes a series for x(i), which he truncates at the term and then 
averages to get an expression for ((x(t))) up to order e^. He then points out 
that this new series is the solution to an ODE, up to order . The stability of 
{{K{t))) is then analyzed in terms of this new ODE. His result for the Mathieu 
equation matches ours up to order (although, he considers the case 7 = 0). 
Our result is a rigorous treatment, applies to higher moments, and we can find 
the solution to any order in e. We stop at in this paper only for convenience. 

If we assume that the Re [A^^^] becomes positive while e is small (so we 
neglect the terms and higher), then we can use ((57|) . ([58]) . and ((59)) to solve 
Re [A(P)] = for p = 1, 2, 3. Then we find that the second moments will become 

unstable before the first moments. If S (^^^AlUq — 7^^ > 5(0), then the third 

moment will become unstable before the second moment. If S ^■\/4wq — 7^^ < 
S{0), then the second moment becomes unstable before the third. 

5.3 Numerical Results 

In this subsection wc discuss the computation of the eigenvalue that determines 
the stability of the Mathieu equation ([T]), with Ao,Ai from (|50p and Fo,Fi 
from ([55)1 . We do not restrict ourselves to small values of e. We carry out these 
calculations by converting the eigenvalue problem to an infinite dimensional 
system of linear equations, and truncating this system after a finite number of 
terms. Our procedure converges rapidly as the number of terms in our expansion 
is increased. 

We limit ourselves to the case of a second-order filter given by ©, with H, 
B, and a given by 

-(T I)' -(J 0)^ -(:). (») 

where /3,ai,a2 G M, 7^ 0, and /J.i,/i2 > 0. The vector of second marginal 
moments m(s), given by ([STjl. satisfies 

Am = Dm + Fom + e(a, s)Fim 

= is^^m-f Mi5.si(sim) - 9^2 ((/?•?! - /^2S2)m) + Fom + e(a,s)Fim. (61) 
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If we multiply (|61l) by and integrate with respect to ds2 , then we get 

Arrij = ^9f^mj + fiids-^isimj) + j/3simj_i + 

- j/i2mj + Fomj + eaisirimj + ea2rimj+i, (62) 

where 

nii(si) = / s^m(si,S2)ds2. 

This is an infinite set of equations for the marginals {mj(si)}. Let (pkisi) = 
-f^fe(v^MiSi)e~^^''S where Hk is the fcth Hermite polynomial. We expand mj in 
the basis (pk as 

mj(si) = ^4v3fe(si). 

fc 

The yjfe are eigenfunctions of the differential operator in the si variable in ([62^ : 
explicitly 

-dg^ipk + /ii9,si(siV3fe) = -kfiiifk fc > 0. 

The Hermite polynomials satisfy the recursion relation Hk+i{y) = 2yHk{y) — 
2kHk-i{y), hence 

siVfc(si) = ^ ^ — {ipk+i{si) + 2fc<pfc„i(si)) fc > 0. 
Thus, simplifies and becomes an equation for 

Ac^^ = (To (fcMi + JM2)I) 4 + ^ {c^Zl + 2(fc + l)c^^j:0 



/Ml 

-^ri (c)-! + 2(fc + l)c^^+i) + ea2ri4+i (63) 



If we consider a finite number of moments rrij for j < N^, and truncate the 
expansion in at fc < Nh-, then we get an approximation to the doubly infinite 
system (|63p . This can be written as a matrix equation 



Lz = Az (64) 

where L is an (NmNhJ) x {NmNhJ) matrix. This eigenvalue problem can be 
solved quickly on a computer. 

Table 1 shows the computed value of A(£) for second moments, which is 
the largest eigenvalue of L in (|64)) . That is, A(e) is the largest eigenvalue for 
the Mathieu equation with filter (in this case the largest eigenvalue is 

real). E2 is the error from a second-order perturbation expansion. That is, 
E2{e) = |Ao + A2e^ — A(e)| with Aq = —7 and A2 is given in equation ([55)1 . i?4 is 
the fourth-order error, i?4(e) = |Ao -I- A2e^ + X4e'^ — A(e)|, where A4 is computed 
by performing the perturbation analysis to order four (the formula for A4 is not 
presented here). The method converges rapidly; the values of A(e) in the table 
were computed for Nm = 7 and Nh = 5. 
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A(e) 


E2 


Ei 


e = 0.01 


-9.89 X 10"^ 


5.74 X 10-« 


2.20 X IQ-^^ 


e = 0.05 


-7.20 X 10-^ 


3.62 X 10-^ 


3.44 X 10~^ 


e = 0.10 


1.65 X 10"^ 


5.96 X 10-4 


2.21 X 10-*^ 



Table 1: Values of the error in computing A(e) (for second moments) for three 
values of e. E2 is the error from the second order expansion, and is the error 
from the fourth-order expansion. Parameter values: /ii — 1.8, /i2 = 0.9,/? = 
1,7 = 0.01, LOO = 0.5, ai = 1, 02 = 0.9, A^™ = 7, iV;, = 5. 



5.4 Alternative Representation of A2 

We present a formula for A2 that involves only Aq and Ai, avoiding construction 
of Fq, Fi. We do not present all of the details because the bookkeeping can be 
quite cumbersome (an interested reader can find the details in [S]), but we 
believe the formula for A2 will be useful for applications. For instance, if one 
wants to compute the perturbation coefficients on a computer, it is easy to build 
an algorithm based on equation (j65|) below, since one only needs to input the 
filter (H,B,a) and the matrices Aq and Ai. 

The equation for the second marginal moments can be written as 

dtM. = Vm + AoM + MAJ + e(a, s) (AiM + MAf ) 

where M is the NxN symmetric matrix with M^j = J^jy XiXjP{s,x.,t)dx.. In 
this case one can solve an eigenvalue problem for the stability where we have 
eigenvalues and eigenmatrices. Looking for solutions of the form M(s,t) = 
e'*'*M(s), yields the eigenvalue problem for M(s) 

AM = VM + AqM + MAj + e(a, s) (AiM + MAJ) . 

The marginal moment tensor M is symmetric {Mjk = Mkj), so we will use a 
basis of symmetric tensors to express M, and in turn reproduce the results of 
The basis that is simplest is given by the eigenmatrices Ejk (and adjoints 
by Fjk with inner product (E,F) = tr(E'^F)) 

E,fe - I {hX + hkhj) , F,k - I (g.gl + gfcgj) , 

where hj are eigenvectors of Ao with eigenvalues aj , and are the normalized 
adjoint eigenvectors of Ao, (gj, h^) = Sjk- The eigenvalues of the Ej^ are sums 
of the a,; AoBjk + Ej-feA^' = (o-j + crfc)Ejfc. 

The analogous result to Lemma [10] is straightforward to show, and following 
the steps in ^we arrive at the following result (note that the eigenvalues i^j of 
Fq from 21 a-nd ^ 35. 21 are sums a-g + <T,n ) ■ 

Theorem 5. Let H and B satisfy the basic conditions (Def. and let {hjjj^^ 
form a complete set. For q, r fixed, if Re [aq + ar — crj — ak + fJ-e] > for each 
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j,k ~ 1, . . . , N and £ = 1, . . . ,n, then the order-two coefficient in the expansion 
A = Aq + A2e^ + . . with Xq = ap + cr,., is given by 



N 




A2 = 8^ 



(65) 



j,k=l 



where 



1 



Cjkhn = - ((5jm(gfc, Aihf) + (5fe™(gj, Aihf) + Sje{gk,Mi^r,i) + 5fc£(gj, Aih^)) , 

and hj are eigenvectors of Aq with eigenvalues aj , and are the normalized 
adjoint eigenvectors of Aq, (gj,hfc) = Sjk- 

6 Conclusions 

We have carried out a perturbation analysis to characterize the moment stabil- 
ity of parametrically forced linear equations, where the forcing is colored noise 
coming out of an Ornstein-Uhlcnbcck process. Our analysis applies to arbitrary 
linear systems, and can in principle be carried out to any order. Our analy- 
sis depends on characterizing the spectrum of the vector Ornstein-Uhlenbeck 
process using ladder operators. Though this spectrum has been characterized 
elsewhere [ITl UHl [22] , we believe the ladder operator approach has been shown 
to be useful in carrying out our perturbation analysis. 
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7 Appendix A: Supplementary Material for §[2] 



In this appendix we give several lemmas used in ^ as well as supplying the 
proofs of several of the lemmas used in that section. 

Lemma 14. The operator V defined in equation ()15|) can be expressed as in 
equation (|2ip . where the djk are the components of the symmetric matrix D, 
given in equation (|23|) . 

Proof. With djk as the components of D given in equation (|23p. we have 

^ 2n-f 1 2?i+l ^ n n ^ 

2 ^ ^ di-jLiLj ~ - ^ ^ {hi-jLiL-j - hi-jL^Lj+n — hj^Li+nLj) — -tr (H) 
i=i j=i 1=1 j=i 



(66) 
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The part of the operator involving the coefficients bij is clearly equal to the 
operator idiv(BV-). To show that the left hand side of equation ([66]) is actu- 
ally T), we need to shows that the terms involving hij are in fact the same as 
~ J27=i K]LiL.j+n = -div (Hs-). We compute 



^ n n ^ 

2 ih,jL,Lj+n + hj,L,+nLj) + ^tr{ii-) 

i=i j=i 

^ n n ^ 

i=i j=i 

^ n 71 ^ n n 

i—l j—1 i—l j — 1 

In the second to last line above, we used the commutator relation from (PO)) . □ 

Proof of Lemma [2] 
Proof. We compute an expression for [D, C] in terms of D and A. 

[2?, /I] ~ ^ ^ '^dijUmiLiLjLjn ^ LjnLiLj) 

= ^ ^ -zdijyjn{Li[Lj , Ljn] + [Li, Lm\Lj) 



2 

2 j,m 



i..j,m i.rii ^ ^ ^i''''^ 

For the equation [X^, = i-iC^ this imphes that we have 



In matrix notation, this is just DAy = /iy, because D = D"^. □ 

This proof holds even if we do not assume that D is symmetric. In that 
case the analysis that follows would be done in terms of the symmetric matrix 
S = ^(D + D-^), instead of D. Thus, it is only for convenience that we use the 
symmetric form of D in ([25)1. 

Proof of Lemma [3] 

Proof. We denote the eigenvalues of H as —fik with Re [/ifc] > for fc = 
1,2, ... ,n. Let u^- be the eigenvectors of H and be the adjoint eigenvec- 
tors 

Hufe = -/ifcUfe, H^Vfe = -p/cVfc (67) 
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■ 










q 






I r 



normalized so that 

(vfc,Uj) = Sjk- 

Recall that H is a real matrix, so complex eigenvalues come in complex conjugate 
pairs. If we write y = (p, q, r)-^ then Ty = yuy becomes 

(68) 

There is a solution with /i = and y° = (0, . . . , 0, 1)-^. If /i 7^ then r = 0, and 
we have two cases. If q = then ([68|) reduces to Hp = /ip. Hence, fi ~ —fik 
and p = Ufe for some k. We will denote this solution as y"*^ = (ufc,O,0)-^. If 
q ^ then we must have H-^q = — ^q, so /i = ^fc and q = Vfc for some k. We 
denote the solution in this case as y*^ = (^(H — /ifeI)~^Bvfe, v^, 0)-^. □ 

Remark 4. T has the eigenvalue 0, with corresponding ladder operator Cq ~ 1. 
This implies that /Zq = 0. However, in this degenerate case, it is convenient for 
notational purposes to define fj,o = — tr (H). We will also write in place of 
—fXk to accommodate negative indices in the proof of Lemma [51 

The Proof of Lemma |4] 

We denote by w^'^'^ the normalized adjoint eigenvectors of T. That is, 
T^w='=''' = ±7Ij,w='='^, (w^'"', y='='^) = Sjk, and (w='='^, y^''') = 0. We begin with a 
preliminary lemma. 

Lemma 15. Let and be the eigenvectors 0/ H as in equations ( [6'7p . Let 
, k ~ 1, . . . ,n, be the eigenvectors of T associated with the eigenvalue ±fJ,k, 
and let •w^'',k = l,...,ri be the normalized adjoint eigenvectors. Then for 
each k = 1, . . . , n, Ay^*^ = w^'^. For fc = 0, 1, . . . , n, Dw='=''' = flkf^^ (using 
liQ = from Remark^. Finally, X]fc=-n^i^2/i ~ ^ij- 

Proof. Note that y^*^ are given explicitly in the proof of Lemma |3] and for fc 7^ 
y-^- = (ufc,O,0)^, y^- = (-(H-/ia)-'Bvfc,Vfc,0)^. (69) 

Wc define 

w'^- = (0, -Ufe, 0)^, w-^ = (vfc, (H - 7IJ)-iBvfc, 0)^ (70) 

and y'' = w'' = (0, . . . , 0, 1)"^. It is straightforward to check that (w*^ , y^'^') = 
Sjk, (w^^y^F*^) = 0, T^wO = 0, and for fc 7^ 0, T^w±''^ = ±7IfeW±^ so 
w^'° are the normalized adjoint eigenvectors. Applying A to the y^'^ in (|69p 
gives Ay='='^ = W^*^ for fc 7^ 0, and hence applying D to Ay*'^ = w^*^ gives 
Dw='='^ = fl^y^^ for fc 7^ 0. With /iq = — tr (H) = /Ig (since H is real), we have 
Dw*^ = Jloy^ ■ (Note that Ay° = 0, so without the convention in Remark |4] we 
would not have Dw*^ = Jloy^.) 

We define the (2?! + 1) x (2?! + 1) matrices Y = [y^", . . . ,y"] and W = 
[w~", . . . , Vif"], then W*Y = l2n+i because (w')'^y^ = Sij for —n < i,j < 
n. But this means YW* = l2ri+i as well, and the components of YW* are 
(YW*)., =EL-n^f2/'- □ 
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We now give the proof of Lemma 01 

Proof of Lemma^ Recall A was defined as having coefficients amp = [Lm, Lp\. 
Writing out [C±j, Ck] in terms of the Lm we have 

2n+l 2n+l 



m,p—l m,p— 1 

= (y±■'")W = (y^^Ay''^ 



Using Ay^*-' = w^*' we Imve (y^^. Ay'"') = (y'''-', w = (y^-?, w^'^). Hence, 
[Cj ,Ck]^ (y+^w-fe) = and , /Ife] = (y-^w-^) = 6,k- □ 

The Proof of Lemma [5] 

Proof of Lemma[B We first consider ^C-kCk = J2l"m=i ^Vni'^UpI'mLp, for 
each k = —n,...,n, using the conventions in Remark 21 For each fc, Dw'^ = 
Mfcy which follows from Lemma [TSl Hence, y"*^ = ^ dmqw'^, so if we 

replace the term y^'' in the above expression for ^C-kCk, and sum over k, we 
get 

n n 2n+l 

^ —C-kCk = Y Y ——dmqwlylLmLp 
k= — n k= — np,m,q=l ^ 

2n+l 



p,m,q=l k= — n 



From Lemmalini ^g^p = V' ^o 



2n+l 2n+l 



^ ^ C,-kCk — ^ ^ '^dmq^qpLmLp — ^ ^ "zdmpLmLp — T> . (71) 



p,m.q—l p.rn—1 



2 '-k'—k — 2 



[C-kCk — 1) by the result of 



LemmajH Combining this with (|7ip and using = ~tr (H) we can write I? as 

V = -Itr (H) + j2[^c^kCk + f - 1)} . 



k=l 



But the eigenvalues of H are —fik, hence tr (H) = — Mfe smd we have 

2^ = J2k=l P-k^-kJ^k- □ 

8 Appendix B: Supplementary Material for 

Proof of Lemma [6l 
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Proof. Suppose x is an eigenvalue of D with eigenfunction (j>, |(/)p(is = 1. If 
wc multiply pT|) by cf), use the definition of T> in ([T5|). integrate over all of space, 
and integrate the term involving B by parts, we get 



^1 



-div (BV0) - 0div (Hs^) ds 
i (V(/), BV0) + 0div (Hs0) ds 



The matrix B is positive semi-definite, so (V(/), BV0) > 0, hence Re [x] < 
Re [— /g„ 0div (Hs0) ds] . But, because H is real, 



2Re 



iiv (Hs0) ds 



iiv (Hs0) + 0div (Hs0) ds. 



If we integrate the first term on the right in this expression by parts, and expand 
the second term we get 



2Re 



(f>div (Hs0) ds 



/ -Wcl) ■ {(j>Us) + 0(0tr (H) + (Hs) • V(/))ds 
: [ |0ptr(H)ds = tr(H) . 

JR" 



Hence, Re [x] < -itr(H). 



□ 



Proof of Lemma [7] 



The proof of Lemma [7] follows almost immediately from a few preliminary 
lemmas. 

Lemma 16. Suppose the eigenvectors of H are complete and the adjoint 
eigenvectors pfe are normalized so (pj,qfc) — Sjk- Let P = [pi, p2, . . . , p^], 
Q [qi, q2, ■ • • , Qn]- We have 



HE 



B 



(72) 



where S ^ = PQ 



Remark 5. Lemmas 1171 and [71 show that, with the appropriate assumptions on 
H and B, the matrix P is invertible and thus there exists a nonsingular matrix 
S = QP^^, so our use of the notation is appropriate. 

Proof. According to equation we have Hp^ + Bq^, = HkPk, and — H-^q = 
A^feQA:- Writing this out in matrix form we get HP + BQ = PM, — H'^Q = QM. 
Here M is the diagonal matrix with fik on the fcth diagonal. Using the second of 
these equations to write M in terms of Q and H, and assuming Q is invertible 
(the eigenvectors of H are complete) we get M = — Q^^H-^Q. Substituting 
this into the first equation we get HP + BQ = — PQ^^H-^Q. If we multiply 
this by on the right and rearrange, we get the result of the lemma. □ 
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We will use the following result for controllable pairs, which follows imme- 
diately from Theorem 2 in [TOj . 

Lemma 17. // B is positive semi- definite, and the eigenvalues of H all have 
real parts less than zero, then the solution to HR + RH"^ = — B is symmetric 
and positive definite provided (H,B) form a controllable pair. 

Lemma [7] follows almost immediately from the previous two lemmas. 

Some lemmas used in the proof of Theorem [2] 

Lemma 18. For any integer m > 0, the operators Ck and C-k satisfy 

[£™+\A] -(m + 1)/:™ (73) 

Proof. For m = 0, this follows immediately from Lemma|31 We can now proceed 
by induction. In particular, if £™j,£fc — = mC™^^, then if we multiply 

both sides of this equation by C-k and use C^^i^CkC-k = C™). {—I + C-uCk) = 
CTit^Hk - we find that - CkCJ^t^ = (m + 1)/:™^, which proves 

the lemma. □ 

Lemma 19. Let H and B satisfy the basic conditions (Def. and let k = 
(fci,fc2,..fc„) be a vector of nonnegative integers. Let 

^^.{s)=C\C%....C\M^), (74) 

then $k(s) is nonzero, and has an eigenvalue of 

n 

Xk = - XI ^Jl^i (75) 

Proof. We begin by showing that C™:^,^o{^) is nonzero for all m > 0. This 
clearly holds for m = by Lemma [51 By induction we can see that if it is 
nonzero for m — 1, then it is non-zero for m. This follows from the fact that 
= to£™^^, and the fact that £fc$o ~ 0. Combining these two facts 
we get ~CkC"2k'^Q{s) = mC^-^^Q{s). This shows that if C"::^^ vanished,then 
-C™^^ would also have to vanish. Since we are assuming this is not the case, it 
follows that £™^,$o(s) does not vanish, and hence by induction does not vanish 
for any m > 0. 

To show that a general function $k(s) does not vanish, we can proceed by a 
different induction proof. In particular, since the operator £_i commutes with 
both £_2 and £2 we see that for any operator Z of the form Z = C^_^ where 
p is a non-negative integer, we have [Z£™2;'^2] = mZC^2^- We can now use 
almost the identical argument as in the last paragraph to show that any function 
of the form Z£"^2^o will be non-zero. We can now carry out this process by 
induction to see that any function of the form $14(8) will be nonzero. 

Once we know that ^kls) is nonzero, it is clear from the ladder operator 
formalism that its eigenvalue must have the form in (|75p . □ 
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There is one subtle point we would like to discuss in our proof of Theorem 
m Our proof relies on the fact that if (f) is an eigenfunction of then either 
Ck(t> ^ 0, or Ck4' gives a new eigenfunction whose eigenvalue has a smaller real 
part. This relies on the assumption that Ck4i remains in the domain of our 
operator. The domain of our operator consists of functions that have moments 
of all orders. Clearly, if this is true of (j), this will be true of £fc0- However, we 
must also make sure that the function Ck(f> has sufficient numbers of derivatives 
to satisfy our differential equation. This is clearly true of the eigenfunctions we 
have found. That is, they clearly have infinitely many derivatives. However, we 
should consider the possibility that there are other eigenfunctions that we have 
not accounted for that are not infinitely differentiable. General theorems on 
elliptic operators rule out such eigenfunctions if B is positive definite. However, 
we have only required that B be positive semi-definite, and that H and B 
form a controllable pair. A heuristic argument that we have found all of the 
eigenfunctions in this less restrictive case is as follows. If we perturb the matrix 
B to make it positive definite, then we know we have all of the eigenfunctions. 
As our perturbation parameter goes to zero, there is nothing unusual happening 
to our spectrum (such as eigenvalues going off to infinity, or clustering about a 
point). Hence, if the eigenfunctions are complete for positive definite B they are 
clearly complete in the less restrictive case where B and H form a controllable 
pair. 

9 Appendix C 

In this appendix, we provide formulas for the asymptotic autocorrelation func- 
tion of the process s{t) and the extended power spectral density (defined in ([75)) ) 
for s{t) as well as for the filter (a, s{t)). In particular, the results of Theorem[S] 
and Corollary [3 are used to express A2 in Theorems [3l HI and [Sj and through- 
out gSl Corollary [7] gives a practical formula for computing the power spectral 
densities of s{t) and (a, s(t)). 

9.1 The Asymptotic Autocorrelation Function 

We begin by proving a lemma concerning the autocorrelation function of s{t) 
as defined in equation s{t) is not a stationary process, but as t — >■ cx) it 
approaches a stationary process, which we refer to as asymptotically stationary. 

Lemma 20. Suppose H and B satisfy the basic conditions (Def. and let 
s{t) be the solution to equation ([2]) with zero initial conditions. As t — > 00 the 
autocorrelation function R(t) = ((s(t)s"^(t + t))) is given by 

R(r) = S^^e"^^ for t > 0, (76) 
and = PQ^^ satisfies equation \7°^ . 
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Proof. We define 



K(t) = e"*Be"^*, Kq = lim / K{t - a)da. (77) 
The solution to equation ^ (with zero initial conditions) is given by 

s{t) = f e^^'~'^^{s)ds 
Jo 



We can write 

ft+T 



s(i)s^(i + r)= / / eH(*--)^(s)|^(r)e""(*+^-'^)drds. 
Jo Jo 

If we take the expected value of both sides of this equation, and use the fact 
that {{^{s)^'^ (r))) = B(5(r — s), we arrive at the equation 



;(s(<)s^(t + T))) = reH(*--)Be""(*-^)e""^da= f K{t ~ a)dae^" 
Jo Jo 



(78) 



When deriving equation (j78p we have assumed that the variable r is equal to 
the variable s at some point when doing the integration. This will only be 
guaranteed if r > 0, and hence this is only valid for r > 0. The expression 
for T < 0, is obtained by using the fact that the autocorrelation function must 
satisfy R(-r) = R^(t). 

Assuming that all of the eigenvalues of H have negative real part, the process 
s{t) will become stationary as i — >■ oo. We take the limit of equation ([75)) as 
t — >■ CX3 to get 

R(t) = KoeH"^ 

where Kq is defined in equation ([77]) . We now show Kq = by showing Kq 
satisfies equation i.e. HKo + KoH^ = -B. 
We have from 

4k(s) = HK(s) + K(s)H'^. 

as 

It follows that 

HKo + KoH^ = - lim / — (K(< - s)) ds. 

t^ocj^ ds 

We can evaluate this integral using the fundamental theorem of calculus. When 
we do this we find that the contribution at s = vanishes in the limit as t — > oo. 
Since K(0) = B, the contribution at s = t is just — B, which completes the proof 
of the lemma. □ 
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9.2 The Extended Power Spectral Density 

The expression for the eigenvalue (with largest real part) of the perturbed op- 
erator V + Tq + e(a, s)ri will be written in terms of the Laplace transform of 
the asymptotic autocorrelation function of the asymptotically stationary filter 
(a, s(<)), which we denote by G. G can be viewed as an extension of the power 
spectral density, and has the advantage that it can be evaluated at points in the 
complex plane, outside of the domain of the power spectral density. 

Definition 7. Let s(t) be an asymptotically stationary stochastic process (i.e. 
stationary in the limit t — > oo) with asymptotic autocorrelation function R(t). 
We define the extended power spectral density of s(i) as 

/•oo 

G(z)= / R(T)e-^^dr. (79) 

With this definition, the scalar filter (a, s{t)) has extended power spectral den- 
sity G{z) = (a, G(z)a). G is indeed an extension of the power spectral density 
S(a;) — /jjj R(r)e~'"'^c?T, because the domain of G contains the set {z G C : 
Re [z] > 0}. In particular, Rc[G(ia;)] = iS(a;), which follows from R"^(t) = 
R(-r). 

Theorem 6. //H andB satisfy the basic conditions (DeJ. then the extended 
power spectral density G(z) for the asymptotically stationary process s[t), de- 
fined in ([2|), is given by 

G(z) = -S-i(H^-zl)~\ (80) 

provided Re [/i; + z] > for I = 1, . . . , n. 

Furthermore, the extended power spectral density G{z) for the asymptotically 
stationary filter (a, s(i)) can be written as 

n „ 

G(z) = (a,G(z)a) = -^-^, (81) 
1=1 ^^^^^ 

where ai, j3i are defined in h8!^) . 

Proof. In Lcmma[2Dl wc showed that the autocorrelation function of the asymp- 
totically stationary process s(i), in the limit t — > oo, is given by R(r) = 
S~^e^^'^ where 

S-i=lim reH(*-^)BeH"(*-^)ds. (82) 
From R(t) = S~^e^ , wc have 

R(i)e"^* dt = -S-i (H^ - zl)^^ , 
assuming that Re [/i; + z] > for Z = 1, . . . , n so that the integral converges. 
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Since a is real, we can use to write a = X]fe=i Q^feVfc = X]fe=i Q^fcVfc- 
Recall, we defined v/ so that H^v/ = — 7I;V/, so we have (H^ — zI)^^Vi = v; 

and 6^^'*""'^^; = e"^'^*~'*-'vi . Using these expressions along with (|86|) . (|82|) . 
and B = B"^, we compute 

G(z) = - hm / y a„azv^eH(*-^) B e^'^*-^' (H^ - ziy^i ds 

n — ft 

= y ^^(v„,Bv;) lim / e-(^'"+'^')(*-^)ds 

□ 

Corollary 7. //H and B satisfy the basic conditions (Def. [^p, then the power 
spectral density S{uj) of the asymptotically stationary filter {a.,s{t)) is given by 
S{ijj) = (a, S{Lu)a) , where S{uj) is the power spectral density of the asymptotically 
stationary process s{t), defined in ([2]), and 

S{uj) = {n^ + iLuiy^ B {n^ ~ iLoiy\ (84) 
Proof. Using the expression for G in equation (|80p we get 
S{uj) = 2Re [G{iiu)] = G{iuj) + G{iuj)* 

= -E-i (H^ - iLoiy^ - (H + iujl)'^ 

= - (H + iiuiy^ ((H + iujl)'E-^ + S-i(H^ - iuji)) {U^ - iLuiy^ 
= - (H + iLoiy^ (HS-i + S-^H^) (H^ - iujiy^ 

= (H. + iLuiy^ B {n^ ~ iujiy^ . 

The asymptotic autocorrelation function for (a, s(t)) = a^s{t) is given by 
((a^s(t)s^(t + r)a)) = (a, R(T)a). Hence S{uj) = (a, S(a;)a). □ 

10 Appendix D: Supplementary Material for 

Proof of Lemma [9] 

Proof. We begin by defining 



afc = (U^a)^,, /3, = -^ ^^^(vfe,Bv„) (85) 



where U = [ui,U2, . . . ,u„]. Recall, {uj} are the eigenvectors of H and {vj} 
are the normalized adjoint vectors. 
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With y='='^ = (p±fc, q±fc, 0)"^, from Lemma [31 we know the ladder operators 
can be written as 

Ck = Pfc • V + Qfc • s, k = -n, . . . , -1, 1, . . . , n, 

with p±k and q±k given exphcitly in the proof of Lemma [31 From these we see 
that for ([37| to be satisfied we must have 

n n n 

^ afcVfc = a, ^ /3fcUfe = ^ aj{H - ^jjI)"^Bvj. (86) 

k=l k=l j=l 

Hence, Vcc ~ a, where V = [vi, V2, . . . , v„]. But V*U — I, so the first expres- 
sion in (|86p is equivalent to the definition of ak in (|55|) . Also, since the {u^} 
are complete, and ((H — fijT)^^)*Vk = —{'Pk +'Pj)^^^k we conclude 

n 

/3fc-(vfc,^a,(H-Ai,I)-iBvj) 

n n 

where the last equality follows from a rearrangement of the sum over j , and the 
fact that the eigenvectors Vj and eigenvalues iij come in conjugate pairs. Thus, 
with ak, Pk defined as in (j85p . the equations in (|86l) arc satisfied, and therefore 
(1371) holds. □ 



Proof of Theorem H 

Proof. From Lemma [T3l we have 

h = -{t^„}_^akT,Ck) = -(Vi, 2^ 2^ ^I'^J^ 

k—1 m—1 J — 1 •' 

J 

= 5^(V'i,ri0^.)(^^,ri</,i)G(;.i - z.,). 

The last equality follows from (|8ip . □ 
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